#include <sys/time.h> 
#include <unistd.h> 
#include <stdlib.h>
#include <stdio.h> 
#include <time.h>   // For srand and rand
#include <string.h> // For memset

int main(int argc, char *argv[]) 
{
    float *a, *b, *c;
    long int i, j, k, size, m;
    struct timeval time1, time2;

    if (argc < 2) { 
        printf("\n\tUsage: %s <Row of square matrix>\n", argv[0]); 
        exit(-1); 
    }

    size = atoi(argv[1]);
    m = size * size;

    // 使用 posix_memalign 分配 16 字节对齐的内存
    if (posix_memalign((void**)&a, 16, m * sizeof(float)) != 0 ||
        posix_memalign((void**)&b, 16, m * sizeof(float)) != 0 ||
        posix_memalign((void**)&c, 16, m * sizeof(float)) != 0) {
        fprintf(stderr, "Memory allocation failed.\n");
        exit(EXIT_FAILURE);
    }

    // 初始化随机种子
    srand(time(NULL));

    // 初始化矩阵 A 和 B
    for (i = 0; i < size; i++) {
        for (j = 0; j < size; j++) {
            a[i*size + j] = (float)(rand() % 1000) / 100.0f;
            b[i*size + j] = (float)(rand() % 1000) / 100.0f;
        }
    }

    // 初始化矩阵 C 为 0
    memset(c, 0, m * sizeof(float));

    // 获取开始时间
    gettimeofday(&time1, NULL);

    // 优化后的矩阵乘法：ikj 循环顺序
    for (i = 0; i < size; i++) {
        for (k = 0; k < size; k++) {
            float a_val = a[i*size + k];
            for (j = 0; j < size; j++) {
                c[i*size + j] += a_val * b[k*size + j];
            }
        }
    }

    // 获取结束时间
    gettimeofday(&time2, NULL);

    // 计算总执行时间
    time2.tv_sec -= time1.tv_sec;
    time2.tv_usec -= time1.tv_usec;
    if (time2.tv_usec < 0L) {
        time2.tv_usec += 1000000L;
        time2.tv_sec--;
    }

    printf("Execution time = %ld.%06ld seconds\n", time2.tv_sec, time2.tv_usec);

    // 释放内存
    free(a);
    free(b);
    free(c);

    return 0;
}